Levy Model for Interstellar Scintillations 
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Observations of radio signals from distant pulsars provide a valuable tool for investigation of 
interstellar turbulence. The time-shapes of the signals are the result of pulse broadening by the 
fluctuating electron density in the interstellar medium. While the scaling of the shapes with the 
signal frequency is well understood, the observed anomalous scaling with respect to the pulsar 
distance has remained a puzzle for more than 30 years. We propose a new model for interstellar 
electron density fluctuations, which explains the observed scaling relations. We suggest that these 
fluctuations obey Levy statistics rather than Gaussian statistics, as assumed in previous treatments 
of interstellar scintillations. 

PACS numbers: 95.30. Qd, 98.38.Am, 98.38. Dq, 95.85.-e 



1. Introduction. Electron density fluctuations in the 
interstellar medium (ISM) cause scintillations of the in- 
tensity of signals arriving from distant pulsars. If the 
medium were completely transparent, the shape of the 
arriving signal would coincide with the shape of the sig- 
nal emitted by the pulsar. However, the observed pulse is 
much broader, and this effect is attributed to the random 
refraction the waves experience while they travel through 
the medium [1-6]. To investigate pulse broadening one 
can assume that the pulsar intrinsic signal is narrow in 
time, Io(t) oc 5(t — to), where Io(t) is the signal inten- 
sity. The observed signal is broad and asymmetric, with 
a sharp rise and a slow decay; see Fig. 1. Observations 
show that broadened shapes of the pulses are similar for 
different pulsars (after proper rescaling), suggesting that 
the density fluctuations statistics along different lines of 
sight are to some extent universal. 

For estimates assume that the pulsar distance is d ~ 
10 kpc, the typical electron density is n ~ 0.03 cm~ 3 , and 
the observational wave frequency is v ~ 500 MHz. Then 
the plasma electron frequency ui pe = (47rne 2 /m e ) 1 / 2 is 
much smaller than v, and density fluctuations change 
the wave phase only slightly. To estimate the time delay 
one can use the approach of geometric optics, where the 
propagating ray is refracted (scattered) by small prisms 
of density inhomogeneities [8-12]. At each scatter event 
occurring at a mean- free path I, the propagation angle 
changes by a small amount, A9 ~ A 2 roAn [see below], 
where A is the wavelength, ro = e 2 /m e c 2 is the classi- 
cal radius of the electron, and An is the density differ- 
ence at characteristic separation I. Using the standard 
assumption that A9 is random and Gaussian, one finds 
that the path direction deviates from a straight line by 
9 ~ A 2 ro Ano(d/l) 1 ' 2 , where Auq is the characteristic 
amplitude of density-difference fluctuations, and the path 
length deviates from the distance d by Ac? ~ d9 2 oc A 4 <i 2 . 
The broadening time can be estimated as T4 ~ Ad/c, 
which gives the standard scaling r<j oc A 4 c? 2 . 

Observations show that the signal width, t^, estimated 
at the half-amplitude level, scales with the wave length 
according to the obtained formula, oc A 4 , while the 



scaling with distance is close to Td oc d , contradicting 
the analytical prediction, as is seen in Fig. 2 [13]. This 
paradox was first discussed by Sutton [1], and although 
the theory of scintillations has been developed for more 
than 30 years, the contradiction has resisted analytical 
understanding [2,3]. 
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FIG. 1. Intensity of a typical observed pulsar signal av- 
eraged over many periods of pulsation. The shown time 
interval spans the pulsar period. The data were taken 
with the Arecibo telescope, at 430 MHz. Courtesy of N. 
D. Ramesh Bhat [7]. 

In this paper we propose that the anomalous scaling 
with the distance is an evidence of non- Gaussian density 
fluctuations in the ISM. We suggest that the probability 
distribution of density gradients has a power-law decay, 
and its second moment is divergent. Such probability 
distributions are common in theories of turbulence, as is 
consistent with the argument that the density statistics 
are governed by turbulent motions in the ISM [14,15]. 
The sum of many angular deviations caused by such 
fluctuations does not have a Gaussian distribution; in- 
stead, the limiting distribution is of the Levy type, and 
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the ray angle performs a Levy flight instead of a con- 
ventional random walk. We present a solvable model 
of scintillations that allows us to unify and extend to a 
non-Gaussian case the standard analytical approaches, 
see, e.g., [8,4]. We then apply this model to Levy den- 
sity statistics, compare it to the observational data, and 
demonstrate that the model naturally produces correct 
scalings of the signals. We report main results here, the 
detailed discussion is presented in [11]. 
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FIG. 2. Pulse temporal broadening as a function of 
the dispersion measure, DM = n(z) dz, which is a 
measure of the distance to the pulsar [16]. Except as 
noted, data were taken from [17]. The solid line has 
slope 4. 

2. Wave equation in a random medium. The Fourier 
amplitude of electric (or magnetic) field, £7 w (r), in the 
isotropic ISM with dielectric permittivity e w obeys the 
wave equation: 



K>,fi(ri,r 2 ) = $ w +n/2(ri)$*_ a/2 (r 2 ), which can be de- 
rived from Eq. (2). Assuming that !l<w,we obtain 



.dV 2k + Akd 2 V 2k-Akd 2 V 2irr 



dz 



4fc 2 9x 2 2 
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AnV, (3) 



lj/c, Ak = Q/c, and An 



where we denoted k 
n(xi, z) - n(x 2 , z). 

Equation (3) is hard to solve without further simplifi- 
cation since n(x, z) is an unknown random function. The 
standard procedure is to assume that the density fluctu- 
ations are Gaussian with a specified correlator in x and 
only short-scale correlations in z, see [4]. Eq. (3) can 
then be averaged over the Gaussian ensemble of density 
fluctuations and over different positions in x-space. How- 
ever, the resulting solution yields a scaling of Td oc \ A d 2 
that contradicts observations, as noted above. 

We propose that the turbulent gas motions in the ISM 
give rise to strongly intermittent and non-Gaussian den- 
sity fluctuations. If the distribution function of An has 
a power-law decay as |An| — > oo and has no second 
moment, then the sum of many independent ray angle 
deviations does not behave as a Gaussian variable (the 
central limit theorem does not hold). Instead, the limit- 
ing distribution, if it exists, is the Levy distribution. A 
random walk whose increments are Levy distributed is 
called a Levy flight. Such processes are common in vari- 
ous random systems, and often replace Brownian motion 
in turbulent systems [19]. 

The Fourier transform (the characteristic function) of 
a symmetric Levy distribution P^An) has the simple 
form, 
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J dAn Pp(An) exp(^An) = cxp ( — C|^| ^) , (4) 



-A - — e w (r) 



E u (t) = 0, 



(1) 



where e a ,(r) = 1 — w 2 e (r)/w 2 , and the electron plasma 
frequency w pe (r) changes slowly on the wave scale A. 
Assuming that the wave propagates in the line-of-sight 
direction, z, we separate the quickly-changing phase of 
the wave from the slowly-changing amplitude, E u (r) = 
exp(iz(jj/c)<f> UJ (z, x), where x is a coordinate perpendicu- 
lar to z. Substituting this into the wave equation (1), we 
derive the equation for the wave amplitude, 



lu d 

2i- — + Aj_ - 47rr n(x, z) 
c dz 



$ w (x,z) = 0, 



(2) 



where Aj_ is a two-dimensional Laplacian in the x plane. 
Following [4-6], introduce the function 7(ri,r 2 ,t) = 
$(ri, t)$*(r 2 , t), whose Fourier transform with respect to 
time is / n (ri,r 2 ) = J dui $ w+ o/ 2 (ri)$* _ n/2 (r 2 ). 

For coinciding coordinates this function is the intensity 
of the radiation whose variation in time we seek. To 
find this function, we may first solve the equation for 



where < (3 < 2, and C is some positive constant. 
Eq. (4) can be taken as the definition of a symmetric Levy 
distribution. The sum of N Levy distributed variables 
scales as An ~ N 1 ^, which becomes diffusion in the 
Gaussian limit (3 = 2. For j3 < 2, the probability distri- 
bution function has algebraic tails, Pp(An) <~ |An|~ 1_/3 
for | An | — ► oo, and its second moment is divergent. We 
thus assume that the random density-gradient fluctua- 
tions are Levy distributed, and are short-scale correlated 
in z. In the next section we first show how Eq. (3) can be 
solved for a general case of non-Gaussian random density 
field. In Sec. 4 we apply our method to Levy distribution. 

3. Batchelor approximation. Propagation as described 
by Eq. (3) cannot be simplified in general, when An is 
not a short-scale correlated Gaussian random variable. 
However, analytical investigation is possible in the im- 
portant case of smooth turbulent fluctuations. This case 
is analogous to the Batchelor limit, in the problem of tur- 
bulent random advection [20] . For that approximation in 
this paper, we neglect all effects other than those of den- 
sity gradients: n(xi) — n(x 2 ) ~ a(z) ■ (xi — x 2 ), where 
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the density gradient c(z) is a random variable with cor- 
relation length I <C d along z. In this approximation, the 
variables separate in Equation (3) and it can be solved 
exactly. We leave analysis of more complicated cases for 
further communcation, and present here results for this 
simple case, which captures the essential physics. 

As a further simplification, consider one-dimensional 
variables xi and x 2 . Since the variables separate, we can 
look for the solution in the factorized form V(x\, xi, z) — 
Ui(xi, z)U~2(x2, z). Then the equation for U\ reads: 



. dUi 

dz 



2k 



Ak d 2 U x 2irr a 

H ; — a(z)xiUi. 



Ak 2 dx\ k 



(5) 



The analogous equation for XJ<i is obtained by changing 
k — > —k. The solution of equation (5) is sought in the 
form U\(xi, z) — A(z) exp [iB(z)x\ + iC(z)xf\ , with the 
initial condition U\(z = 0) = S(x\), if the refracting 
medium extends all the way up to the pulsar. Substi- 
tuting this ansatz into (5), we find 



A(z) — —= exp 



-i(2k-Ak) 
4fc2 



B 2 {z')dz' 



Z 

B(z) = J a(z')z'dz', C(z) 



(2k - Ak)z 



(6) 



(7) 



Note that this solution describes the path of a single 
ray through a sequence of density gradients <j(z). Ef- 
fects of multiple rays can be found from superposition. 
The intensity of received radiation can be calculated 
from the Fourier transform I u (z,t) = J^° 00 V j j i q(x — 

0, z) exp(— iSlt)dSl/\/2n. In this Fourier transform, indi- 
vidual ray paths will yield contributions with phase pro- 
portional to f2 = Akc, with coefficient equal to the travel 
time for that path. Cross-terms describing interference of 
paths yield contributions that oscillate rapidly with fre- 
quency and average to zero, see e.g. [12]. The intensity, 
averaged over an ensemble of statistically independent 
rays, is then given by the average over individual travel 
times or equivalently over different realizations of <j(z). 
This leads to 



B 2 (z')dz' ), 



(8) 



where the angular brackets denote the statistical average. 
Note that B(z) is proportional to the deflection angle 9 
of the ray. Formula (8) gives the shape of the signal 
observed at the Earth; if the scattering medium were ab- 
sent, the signal would be undistortcd, I(t) oc 5(t — d/c). 

Lee & Jokipii investigated Eq. (3) for short-scale cor- 
related Gaussian density fluctuations [4]. For averaging 
over Gaussian <r(z), our solution (8) reproduces those 
obtained by Williamson [8] with a phenomenological ap- 
proach. Thus, Williamson's solutions are applicable un- 
der the assumption of smooth Gaussian density fluctu- 
ations, and when one keeps only the linear term in the 



expansion of An. In the next section we apply our ap- 
proach to the Levy distributed density fluctuations. 

4. Scintillations as Levy flights through the interstel- 
lar medium. The averaging in formula (8) can be per- 
formed for the Levy distributed short-scale correlated 
density gradients a(z). To do this, we represent inte- 
grals is (7), and (8) in the discretized forms, i.e., we 
assume that d = nl, z' = ml, and z" — si, where 
I is the correlation length of density fluctuations, and 
change J* f(z')dz' — > J2m=i f(^ m )^ f° r an arbitrary 
function f(z). The right hand side of (8) is the probabil- 
ity distribution of the propagation time delay, r = t—d/c. 
For a continuous medium, this time delay is given by 



r 2/3\4 ™ 

T ~ 8tt 2 c ^ 



1 m 

m z — ' 



sa s 



(9) 



The solution (8) can now be calculated numerically as 
the probability distributions of this variable, under the 
assumption that a s are distributed independently, iden- 
tically, and according to the Levy law (4). We how- 
ever need to specify the parameter (3 in the Levy for- 
mula. We do this by comparing our model with observa- 
tions. The observed scaling of pulse broadening is close 
to Td oc A 4 rf 4 , while our model gives r cx 
as is seen from the scaling for sums of Levy-distributed 
variables, ^ m a s <~ to 1 /^ 3 , following Eq. (4). Thus, we 
obtain, (3 w 2/3. 

Note that standard Gaussian models of density distri- 
bution were not able to satisfy both observational scal- 
ings, A 4 and c? 4 , simultaneously. Various studies of non- 
smooth density fluctuations within these models have not 
reproduced this scaling cither [10,11,21]. Our model re- 
produces the anomalous ensealing naturally. Moreover, 
it predicts that the probability distribution function of 
electron density gradients in the interstellar turbulence 



decays as P(a) <~ |a| 



-1-/3 



led 5 / 3 . Power-law distribu- 



tions P(<t) with (3 < 2 are indeed observed in numerical 
simulations of compressible turbulence [22], however, no 
one has yet derived them from first principles. To date 
theories of scintillations have exploited only second-order 
correlators of the density fluctuations, while in our ap- 
proach these correlators do not exist (or do not matter) 
and one must work with the whole probability distribu- 
tion function. 

Although our goal was to explain the scalings of the sig- 
nals, it is interesting to see to what extent we can predict 
their shapes. The delay time is proportional to the square 
of the typical deflection angle of the ray trajectory, r oc 
8 2 , where 9 has the Levy distribution Pp(9). Therefore, 
the distribution of arrival times is I(t) oc P ( 3(t 1 / 2 )t -1 / 2 , 
with the asymptotic form 7(r) oc t" 1 "' 3 / 2 as r — > oo. 
Fig. (3) shows the distribution of r from numerical cal- 
culation, with a power-law decay at long times, as ex- 
pected. Because the observed shape of the scattered 
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pulse is directly related to the probability distribution 
of gradients in electron density in the ISM, observational 
data offer the possibility of characterizing interstellar tur- 
bulence [7]. 
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FIG. 3. Pulse-broadening functions for the model of 
linear density fluctuations obeying the Levy statistics 
with (3 = 2/3. The distribution of the time delay (9) is 
found numerically using 10 6 rays. The insert shows the 
large-time asymptotics of the curve in the log-log scale. 

The curve in Fig. (3) closely resembles the observed sig- 
nals, although the presented analytical shapes are the re- 
sult of averaging over an ensemble of non-interfering rays, 
corresponding to an observational average over an infinite 
amount of time. In practice, the averaging time is finite, 
and the long tail of the distribution, dominated by rare 
events, may not have converged. We also ignore instru- 
mental response [7]. Moreover, a non-analytic density 
field is more natural for a turbulent cascade [15,14,23]. 
Also, the small-scale density fluctuations that produce 
scintillations should be collisionless, and elongated along 
the local magnetic field [14]; the scattering is likely to 
be highly anisotropic, and locally nearly 1-dimensional. 
We will consider these effects in future work. Interest- 
ingly, some scattered pulsars show power-law declines 
at long times, such as that in Fig. 1. Scintillation of 
nearby pulsars also shows evidence for weak large-angle 
scattering [24]. Some interferometric studies suggest a 
"halo" surrounding the source at large scattering angles 
and excess scattering at small angles relative to a Gaus- 
sian [25,26], as might be expected for a Levy distribution 
of scattering angles. Intrinsic source structure, and the 
relatively short observatonal averages, may complicate 
this interpretation. 

Finally, we wish to comment on the original expla- 
nation of the anomalous ensealing by Sutton [1]. Sut- 
ton suggested that encounters with much more-strongly- 
scattering HII regions become more probable on longer 



lines of sight. This however requires a perhaps surpris- 
ingly close coordination of DM (over 1.7 orders of mag- 
nitude) with Td (over 8 orders of magnitude). Sutton's 
proposal assumes essentially non-stationary statistics for 
the density distribution along z. Our proposal also in- 
vokes rare, large events, but in a statistically stationary 
way. 

To summarize, we propose that the observed anoma- 
lously strong time-broadening of pulsar signals is evi- 
dence for non-Gaussian distribution of electron density 
gradients in the ISM. We argue that this distribution is 
of the Levy type, in accord with the turbulent origin of 
density fluctuations, and we present a simple model that 
explains the observational scalings of pulsar signals. 

The work of SB was supported by the ASCI Flash Cen- 
ter at the University of Chicago, under DoE subcontract 
B523820. 
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